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Solving Einstein's constraint equations for the construction of black hole initial data requires 
handling the black hole singularity. Typically, this is done either with the excision method, in 
which the black hole interior is excised from the numerical grid, or with the puncture method, in 
which the singular part of the conformal factor is expressed in terms of an analytical background 
solution, and the Hamiltonian constraint is then solved for a correction to the background solution 
that, usually, is assumed to be regular everywhere. We discuss an alternative approach in which the 
Hamiltonian constraint is solved for an inverse power of the conformal factor. This new function 
remains finite everywhere, so that this approach requires neither excision nor a split into background 
and correction. In particular, this method can be used without modification even when the correction 
to the conformal factor is singular itself. We demonstrate this feature for rotating black holes in the 
trumpet topology. 
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I. INTRODUCTION 

Constructing initial data in general relativity requires 
solving the constraint equations of Einstein's field equa- 
tions (see, e.g., [2 II]). Under the assumption of confor- 
mal flatness and maximal slicing, solutions to the mo- 
mentum constraint describing boosted or spinning back 
holes can be expressed analytically in terms of Bowen- 
York solutions [3]. These solutions can then be inserted 
into the Hamiltonian constraint, which, in general, has 
to be solved numerically for the conformal factor. 

Any numerical method employed to solve the Hamilto- 
nian constraint for black hole data has to accommodate 
the presence of black hole singularities. One approach is 
the excision method, in which the black hole interior is 
excised from the numerical grid, and suitable boundary 
conditions are imposed on the black hole horizon (see, 
e.g., IIHS])- An alternative is the puncture method, in 
which the singular parts of the solution are expressed 
in terms of an analytical background solution, and the 
Hamiltonian constraint is solved for a regular correction 
to the background solution [7| |H] • 

Here we discuss an alternative approach that requires 
neither excision nor a decomposition into background 
and correction (even though the latter is probably de- 
sirable in terms of numerical accuracy). Specifically, we 
consider solving the Hamiltonian constraint for an inverse 
power of the conformal factor. This approach, which is 
similar to an approach that has become extremely suc- 
cessful in solving Einstein's evolution equations (e.g. [S]), 
appears to be promising in the context of Einstein's con- 
straint equations as well. The new function representing 
the conformal factor remains finite everywhere and can 
be solved for directly. We present numerical examples 
and compare with both analytical and independent nu- 
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merical results. An important advantage of this approach 
over the puncture method is that it can be used without 
modification even when, in the puncture method, the cor- 
rection diverges as fast as the background solution itself. 
We demonstrate this feature for rotating black holes in 
the trumpet topology TU', TT] . We expect that this prop- 
erty may be important for applications that relax the 
assumption of conformal flatness, since a non-vanishing 
deviation from conformal flatness may lead, in the con- 
text of the puncture method, to singular corrections to 
the analytic background terms. 

II. BASIC EQUATIONS 

In vacuum, and under the assumption of maximal slic- 
ing and conformal flatness, the Hamiltonian constraint 
reduces to 

D'^ = -^^l;-'A,A'^, (1) 

where ip is the conformal factor and A'^^ — ^p^'^A^^ is the 
conformally rescaled, trace-free part of the extrinsic cur- 
vature. Also, Di is the covariant derivative operator with 
respect to the conformally related metric jij = ip'^jij, 
where jij is the physical spatial metric. Under the as- 
sumption of conformal flatness, ^ij = rjij, where r]ij 
is the fiat metric any coordinate system, the operator 
= DiDj reduces to the flat Laplace operator. 
Now consider a new function 

= (2) 

where n is a constant that we will later choose to be a 
positive integer. Since the conformal factor %j) typically 
diverges at the black hole singularity, this choice makes Q. 
go to zero and remain finite there. The Laplace operator 
acting on U, then satisfies 

D^n = ( 1 + - ) n-^fw.nDjfi - nn^+^/'^D^i;. (3) 
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Inserting ([T]) we can now express the Hamiltonian con- 
straint in terms of fl, 



1\ D'nD,n n 1 , o/„ - -,„■ 



Evidently we recover the Hamiltonian constraint in its 
original form ([T]) for n = —1. 

Different choices for the power n can be considered. 
One possibility would be to choose n in such a way that 
that the first term on the right-hand side of Q remains 
finite at the black hole singularity, where Q vanishes. 
Assuming that ip diverges with r"™ at the singularity, 
where r is the isotropic radius, J7 scales with r"™, and 
the first term on the right-hand side of Q with 7-"™-2^ 
This term remains finite as r — >■ if rt > 2/r7i. For the 
minimum value, n — 2 /to, we always have 17 oc close 
to the singularity. 

The scaling of -0 close to the singularity depends on the 
slicing of the black hole. For so-called wormhole data, for 
example, ip diverges with r^^, so that m = 1. Here we 
will focus on trumpet data, for which ip cx r^^/^ close 
to the singularity (see p!^HT5] and equation Q below). 
The above argument would then suggests that we should 
choose n > 4. Despite these considerations, n — A may 
not be the best choice. As we will discuss in more detail 
below, we found better results for n = 2, even though 
some terms in Q diverge with 1/r for this choice. 

Before proceeding we linearize the Hamiltonian con- 
straint Q as follows. Denoting the approximate solution 
after I iteration steps with fi', we search for a correction 
50. ^ il' so that r2'+^ = fi' -I- Jil solves the equation. 
Denoting the residual of equation Q for 0} with 7?.' , 



7^' = D^n^ - ( 1 



\{0'Y+^/^A\ (5) 



where we have abbreviated — AijA^^, the linearized 
Hamiltonian constraint Q becomes 

D^{Sn) - f 1 + - ) ^fWiQ^D^Sn + 



n J 9} 

+ - 
n 

= -n 



A^Un 



(6) 



Equation (|6| is a linear equation that can be solved it- 
eratively for SO until the norm of the residual 7?.' has 
dropped below a desired tolerance. 



III. NUMERICAL EXAMPLES 

A. Schwarzschild 

As a first test of this scheme we solve the Hamilto- 
nian constraint Q in spherical symmetry to construct 
the Schwarzschild solution. 
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FIG. 1: Numerical solutions Q for a Schwarzschild black hole, 
for n = 2 and n = 4. Here we imposed the outer boundary at 
-Rout ~ IGM and used lOOA'-l-l gridpoints. The upper panels 
shows the solution for different values of N. As expected, Q 



scales with 



V2 



at the center. The lower panel shows the 



rescaled errors N'^AQ, demonstrating second-order conver- 
gence to the analytical solution, even in the neighborhood of 
the singularity. 



We adopt a finite-difference method and use a uniform 
vertex-centered grid so that the first grid point is at r = 
0. We set i7 = at r = 0, and set to its analytical 
value ■(/'o'" 3't the outer boundary i?out of the grid. Finite- 
differencing the operator on the left-hand side of equation 
([g]) results in a tridiagonal matrix that can be solved with 
elementary methods. 

In order to construct maximally sliced trumpet-data, 
we adopt 



^0 — 



3\/3M2 



(7) 



where rt* = x^/r is the spatial normal vector pointing 
away from the center of the black hole at r = 0, and 
M is the total mass-energy of the black hole [THj. The 
analytical solution for ipo can also be given analytically, 
albeit only in parametric form (see [II]). In the neigh- 
borhood of the singularity, ipo is given by 



/3A/ 



1/2 



(8) 



to leading order in r. Throughout this paper, Af , denotes 
the mass of the background Schwarzschild solution. 

In Figure [T] we show numerical results as a function of 
the radius r for two choices n = 2 and n = 4. The upper 
panel shows the solutions Q for different grid resolutions; 
as the resolution increases, the numerical solutions ap- 
proach the analytical solution fip = V'o'"- expected. 
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FIG. 2: Numerical solutions w for boosted black holes with 
momentum P" /M = 1.0. The lines in the large plot show 
results for w, obtained for n = 2, on three different numerical 
grids with increasing grid resolution and more distant outer 
boundaries. For A'^ = 65, the outer boundaries were imposed 
at Xout = 16M, for N = 129 at Xout = 24M, and for iV = 257 
at Xout = 32M. Also included, as crosses, are the results as 
computed from the Hamiltonian constraint in its original from 
([l]) with A'^ = 257. The insert shows higher-resolution results 
(obtained with N = 257 and Xout ~ 2M) of the region around 
the black hole center. This graph also includes the leading- 
order analytical result ( 13 1 as a solid line. 



scales with r"/^ close to the center. The lower panel 
shows the errors Ail = U — ^Iq, rescaled with the square 
of the grid-spacing. These results demonstrate that the 
scheme is second-order accurate even in the vicinity of 
the black hole singularity. They also demonstrate that 
this method can be adopted without decomposing the 
solution into a background and a correction term, even 
though we will use such a decomposition in the next sec- 
tions. 

The choice n = 4 has the appealing features that the 
first term on the right hand side of Q remains finite, as 
we discussed above, and that Q is smooth at the center. 
However, Fig.[l]shows that for n = 2 the numerical errors 
are smaller. We have also found that for n = 4 our 
iterative scheme failed to converge for rapidly spinning 
black holes (see Section IIIC below). We will therefore 
focus on ri = 2 for the remainder of this paper. 



B. Boosted black holes 



We construct boosted trumpet black holes by adding to 
the Schwarzschild extrinsic curvature ([t]) a Bowen-York 
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{f^ - n'n^)nkP^) (9) 



(since the momentum constraint is linear, the sum of two 
solutions is still a solution). In the interest of numerical 
accuracy, we decompose f2 as 



(10) 



where fio is the analytical Schwarzschild solution, and w 
a correction. This way the largest terms in the Hamil- 
tonian constraint can be computed from the analytical 
solution, and only the correction w needs to be treated 
with finite-differencing. The iteration for w still works 
in the same way as before if we replace in (joj) with 



w 



and 5VL with 5w. 



We now solve Q, using the iteration ([6]), on a three- 
dimensional, uniform cartesian grid of gridpoints, 
with the help of both Cactus and PETSc software. We 
again adopt a vertex-centered grid, fix w = at the 
center, and impose a l/r fall-off condition at the outer 
boundaries on a square x^nt = ^out = ^out = i-'^out- 

Boosted trumpet black holes have previously been con- 
structed with the puncture method by solving the Hamil- 
tonian constraint in its original form ([T]) (see jlOl [11] , see 
also [16^ for a calculation using the excision method). In 
that approach, the conformal factor ijj is decomposed as 
ip = tpo + u. Given a solution u, we can compute the 



corresponding w from (10), 



w ^ fl — flo ~ ("00 + u) 



(11) 



For boosted black holes, we can therefore compare the 
results from the new method discussed here with inde- 
pendent numerical results. Moreover, as shown in [lOj . 
regular solutions for u for boosted trumpet black holes in 
the neighborhood of the puncture are given, to leading 
order in r, by 
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(12) 



where cos 9 = n^P*. Inserting this, together with (18 



into (11) we find 



2^/2 p 



Wp 



3(n+3)/2 ^/ VA/ 



■i/2+l 



COS 9 



(13) 



to lowest order in r. 

In Fig. [2] we show numerical results for a black hole 
boosted with a momentum P^/M = 1.0. The graph 
shows that, as both the resolution and the distance to 
the outer boundaries are increased, the results for w ap- 
proach those computed with the puncture method from 
u. Note that w and u feature different asymptotic be- 
havior as r — )■ oo. By imposing a l/r fall-off on both 
functions at a finite Xout, we suppress different higher- 
order terms. Therefore, the results from the two ap- 
proaches only approach each other as both the numerical 
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FIG. 3: Solutions w for a black hole spinning with angular 
momentum 3^ jM"^ = 0.1, using n = 2. We show numer- 
ical results (dots) in the vicinity of th e sin gularity, together 
with the leading-order analytical result ( 16 1 (solid lines), both 
along the direction of the spin (cos d = 1) and a direction or- 
thogonal to the spin (cosS = 0). The numerical results were 
obtained with iV = 257 and Xout = 2M. 




FIG. 4: Numerical results for w for a black hole spinning with 
angular momentum J'' /M'^ — 1.0 in the x-z plane, obtained 
with N = 129, Xout = 4A/ and n = 2. 



resolution and the location of the outer boundaries are 
increased. The inset in Fig. [2] shows higher-resolution 
results for the region around the singularity, and demon- 
strates that the numerical results, for both w and u, ap- 
proach the analytical result ( 13 ) as r — )• 



m 



C. Spinning black holes 

Spinning black holes can be constructed by adding to 
the Schwarzschild extrinsic curvature ([?]) a Bowen-York 
solution representing a black hole with spin angular mo- 
mentum J', 



JkTll. 



(14) 



As a consequence of A^g diverging with r ^ (rather than 



now diverge with r"-*^/^, i.e. they diverge with the same 
power of r as the background ■00 itself. In the neighbor- 
hood of the singularity, corrections u, to leading order in 
both r/M and J /M'^ , are given by 



us 



12 V M2 y V 3r / 



1/2 



(15) 



where cos 9 = j* (see pH]). 

Given that the corrections u are not regular, the punc- 
ture method in its original form breaks down. The equa- 
tions can still be solved if the singular behavior of u is 
scaled out, but this leads to rather complicated expres- 
sions (see |11|). The method proposed here, however, can 
still be used without any modi ficat ion (other tha n u sing 
instead of A^). Inserting (15) and ^ into (11) we 
now have 
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for the boosted solutions 



corrections u to ■0o 



to leading order in both r/M and J/M'^. Not surpris- 
ingly, this solution scales with r"/^, just like the back- 
ground term ilp- 

In Fig. [3] we show high-resolution results, for n = 2, 
in the vicinity of the center for a black hole with angu- 
lar momentum J^/M^ = 0.1. For this sufficiently small 
value of J, the numerical results (crosses and circles) ap- 
proach the center as predicted by the analytical result 
(solid lines, see eq. (16)). Note that w remains finite 
everywhere, even if, for n — 2, derivatives are discontin- 
uous at r = 0. However, since we can set the solution 
to zero there, and never need to evaluate any derivatives 
at the center, this does not affect the numerical scheme. 
In Fig. |4] we also show a surface graph of w for a larger 
value of the angular momentum [J'^ /M"^ = 1.0). 



IV. DISCUSSION 

We discuss an approach to solving the Hamiltonian 
constraint that requires neither excision nor a decompo- 
sition into background and correction terms for the con- 
struction of black hole initial data. Specifically, we solve 
the Hamiltonian constraint for an inverse power of the 
conformal factor. The resulting function then remains 
regular everywhere, and vanishes at the location of the 
black hole singularities. We present numerical examples 
and compare with both analytical and independent nu- 
merical results. An important advantage of this method 
is that it can handle cases for which, in the puncture 
method, corrections diverge at the singularity as fast as 
the background term itself. We demonstrate this feature 
for spinning trumpet black holes, and expect that this 
property may be important for applications in which the 
assumption of conformal flatness is relaxed pU] . 

We also experiment with different powers n in the 
rescaling Q = -0"" . As we discuss in Section the 
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choice n = 4 leads to some appealing properties of the 
equation and the solutions; we nevertheless found n = 2 
more suitable for our iteration scheme. 

In the puncture method, corrections u that diverge at 
the singularity, if they exist, are suppressed by the as- 
sumption of regularity there [21j . In the approach dis- 
cussed here, such solutions would not be suppressed au- 
tomatically, since they also satisfy w = at the center. 
We therefore suspect that a second branch of solutions, 
corresponding to singular corrections u, might exist. An- 
alyzing the uniqueness of solutions, and the properties of 



any other branches of solutions, might make an interest- 
ing subject for future investigations. 
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